R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

library(spdep)
## Loading required package: sp
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Loading required package: sf
## Linking to GEOS 3.5.1, GDAL 2.2.2, PROJ 4.9.2
library(maptools)
## Checking rgeos availability: TRUE
library(leaflet)

chi.poly <- readShapePoly('./shapefiles/PovertyData.shp')
## Warning: readShapePoly is deprecated; use rgdal::readOGR or sf::st_read
class(chi.poly)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
str(slot(chi.poly,"data"))
## 'data.frame':    88 obs. of  27 variables:
##  $ OBJECTID_1: int  1 2 3 4 5 6 7 8 9 10 ...
##  $ COUNTY_CD : Factor w/ 88 levels "ADA","ALL","ASD",..: 36 37 38 39 25 40 41 42 61 63 ...
##  $ COUNTY_SEA: Factor w/ 88 levels "AKRON","ASHLAND",..: 33 42 51 58 20 35 70 53 10 61 ...
##  $ ODOT_DISTR: int  9 10 11 3 6 9 11 5 10 1 ...
##  $ FIPS_COUNT: Factor w/ 88 levels "39001","39003",..: 36 37 38 39 25 40 41 42 61 63 ...
##  $ POP_2010  : int  43589 29380 42366 59626 1163414 33225 69709 60921 14645 19614 ...
##  $ POP_2000  : int  40875 28241 38943 59487 1068978 32641 73894 54500 14058 20293 ...
##  $ POP_1990  : int  35728 25533 32849 56240 961437 30230 80298 47473 11336 20488 ...
##  $ STATE_PLAN: Factor w/ 2 levels "N","S": 2 2 1 1 2 2 1 1 2 1 ...
##  $ ELEVATION_: int  1340 1220 1380 1200 1130 1040 1420 1420 1340 780 ...
##  $ ELEVATION1: int  700 640 800 600 670 600 640 840 660 650 ...
##  $ LAT_NORTH_: num  39.4 39.7 40.7 41.3 40.1 ...
##  $ LAT_SOUTH_: num  39 39.4 40.4 41 39.8 ...
##  $ LONG_EAST_: num  -83.3 -82.2 -81.6 -82.3 -82.8 ...
##  $ LONG_WEST_: num  -83.9 -82.7 -82.2 -82.8 -83.3 ...
##  $ AREA_SQMI : num  558 424 424 496 544 ...
##  $ AREA_ID   : int  53 54 55 56 143 57 58 59 75 76 ...
##  $ SHAPE_STAr: num  2.41e+09 1.84e+09 1.91e+09 2.27e+09 2.40e+09 ...
##  $ SHAPE_STLe: num  213158 215452 193045 197349 211754 ...
##  $ NoOfUni   : int  0 0 0 0 13 0 1 2 0 0 ...
##  $ Shape_Leng: num  1.72 1.75 1.59 1.58 1.7 ...
##  $ Shape_Area: num  0.151 0.115 0.117 0.138 0.149 ...
##  $ FIPS      : int  39071 39073 39075 39077 39049 39079 39081 39083 39121 39125 ...
##  $ RUC_code  : int  6 1 7 4 1 7 3 4 7 6 ...
##  $ Pov_low_bo: num  13.9 11.5 7 12.9 15.2 14.3 14.7 8.5 12.7 8.1 ...
##  $ Pov_Upper_: num  19.7 17.5 11 16.7 16.8 21.5 20.5 12.9 19.7 12.3 ...
##  $ Pov_Percen: num  16.8 14.5 9 14.8 16 17.9 17.6 10.7 16.2 10.2 ...
##  - attr(*, "data_types")= chr  "N" "C" "C" "N" ...
str(slot(chi.poly,"bbox"))
##  num [1:2, 1:2] -84.8 38.4 -80.5 42
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:2] "x" "y"
##   ..$ : chr [1:2] "min" "max"
str(slot(chi.poly,"proj4string"))
## Formal class 'CRS' [package "sp"] with 1 slot
##   ..@ projargs: chr NA
str(slot(chi.poly,"data"))
## 'data.frame':    88 obs. of  27 variables:
##  $ OBJECTID_1: int  1 2 3 4 5 6 7 8 9 10 ...
##  $ COUNTY_CD : Factor w/ 88 levels "ADA","ALL","ASD",..: 36 37 38 39 25 40 41 42 61 63 ...
##  $ COUNTY_SEA: Factor w/ 88 levels "AKRON","ASHLAND",..: 33 42 51 58 20 35 70 53 10 61 ...
##  $ ODOT_DISTR: int  9 10 11 3 6 9 11 5 10 1 ...
##  $ FIPS_COUNT: Factor w/ 88 levels "39001","39003",..: 36 37 38 39 25 40 41 42 61 63 ...
##  $ POP_2010  : int  43589 29380 42366 59626 1163414 33225 69709 60921 14645 19614 ...
##  $ POP_2000  : int  40875 28241 38943 59487 1068978 32641 73894 54500 14058 20293 ...
##  $ POP_1990  : int  35728 25533 32849 56240 961437 30230 80298 47473 11336 20488 ...
##  $ STATE_PLAN: Factor w/ 2 levels "N","S": 2 2 1 1 2 2 1 1 2 1 ...
##  $ ELEVATION_: int  1340 1220 1380 1200 1130 1040 1420 1420 1340 780 ...
##  $ ELEVATION1: int  700 640 800 600 670 600 640 840 660 650 ...
##  $ LAT_NORTH_: num  39.4 39.7 40.7 41.3 40.1 ...
##  $ LAT_SOUTH_: num  39 39.4 40.4 41 39.8 ...
##  $ LONG_EAST_: num  -83.3 -82.2 -81.6 -82.3 -82.8 ...
##  $ LONG_WEST_: num  -83.9 -82.7 -82.2 -82.8 -83.3 ...
##  $ AREA_SQMI : num  558 424 424 496 544 ...
##  $ AREA_ID   : int  53 54 55 56 143 57 58 59 75 76 ...
##  $ SHAPE_STAr: num  2.41e+09 1.84e+09 1.91e+09 2.27e+09 2.40e+09 ...
##  $ SHAPE_STLe: num  213158 215452 193045 197349 211754 ...
##  $ NoOfUni   : int  0 0 0 0 13 0 1 2 0 0 ...
##  $ Shape_Leng: num  1.72 1.75 1.59 1.58 1.7 ...
##  $ Shape_Area: num  0.151 0.115 0.117 0.138 0.149 ...
##  $ FIPS      : int  39071 39073 39075 39077 39049 39079 39081 39083 39121 39125 ...
##  $ RUC_code  : int  6 1 7 4 1 7 3 4 7 6 ...
##  $ Pov_low_bo: num  13.9 11.5 7 12.9 15.2 14.3 14.7 8.5 12.7 8.1 ...
##  $ Pov_Upper_: num  19.7 17.5 11 16.7 16.8 21.5 20.5 12.9 19.7 12.3 ...
##  $ Pov_Percen: num  16.8 14.5 9 14.8 16 17.9 17.6 10.7 16.2 10.2 ...
##  - attr(*, "data_types")= chr  "N" "C" "C" "N" ...
summary(chi.poly@data$NoOfUni)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.500   1.295   1.000  14.000
plot(chi.poly)

leaflet(chi.poly) %>%
  addPolygons(stroke = TRUE, fillOpacity = 0.5, smoothFactor = 0.5) %>%
  addTiles() #adds a map tile, the default is OpenStreetMap
chi.poly@data$Pov_low_bo
##  [1] 13.9 11.5  7.0 12.9 15.2 14.3 14.7  8.5 12.7  8.1 16.0 10.2 15.3  6.3 14.6
## [16]  7.5  9.6  5.3  8.1 14.9  5.0 13.3 11.2  9.1 13.6 16.9 18.3 25.7 11.5  7.3
## [31]  8.3 13.1 14.5 11.8 13.5  9.0 10.8 13.5  7.5  7.4 16.5  7.3  9.1 12.1 16.6
## [46]  7.5 16.8 14.0  5.0 16.4 11.8 17.0  8.8 12.4  7.3  7.6 12.9  6.7 15.2  5.4
## [61]  8.3 15.2  7.9 13.0  6.4  4.0 10.7  9.8 12.7 13.2 12.3 13.2 11.1 17.3  7.5
## [76]  8.8  9.4  7.6  4.1 10.3  7.2 12.9 10.5 15.7  3.8  9.7  9.6  9.4
require(RColorBrewer)
## Loading required package: RColorBrewer
qpal<-colorQuantile("OrRd", chi.poly@data$Pov_low_bo, n=4) 
leaflet(chi.poly) %>%
  addPolygons(stroke = TRUE, fillOpacity = .8, smoothFactor = 0.2, color = ~qpal(Pov_low_bo)
  ) %>%
  addTiles()
chi.ols<-lm(  NoOfUni~Pov_low_bo, data=chi.poly@data)
summary(chi.ols)
## 
## Call:
## lm(formula = NoOfUni ~ Pov_low_bo, data = chi.poly@data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2727 -1.2258 -0.6412  0.1768 11.8641 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.20538    0.81222  -0.253   0.8010  
## Pov_low_bo   0.13534    0.06905   1.960   0.0532 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.542 on 86 degrees of freedom
## Multiple R-squared:  0.04276,    Adjusted R-squared:  0.03163 
## F-statistic: 3.842 on 1 and 86 DF,  p-value: 0.05323
list.queen<-poly2nb(chi.poly, queen=TRUE)

list.queen
## Neighbour list object:
## Number of regions: 88 
## Number of nonzero links: 460 
## Percentage nonzero weights: 5.940083 
## Average number of links: 5.227273
W<-nb2listw(list.queen, style="W", zero.policy=TRUE)

W
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 88 
## Number of nonzero links: 460 
## Percentage nonzero weights: 5.940083 
## Average number of links: 5.227273 
## 
## Weights style: W 
## Weights constants summary:
##    n   nn S0       S1       S2
## W 88 7744 88 35.22186 356.7656
plot(W,coordinates(chi.poly))

coords<-coordinates(chi.poly)
W_dist<-dnearneigh(coords,0,1,longlat = FALSE)


moran.lm<-lm.morantest(chi.ols, W, alternative="two.sided")

print(moran.lm)
## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = NoOfUni ~ Pov_low_bo, data = chi.poly@data)
## weights: W
## 
## Moran I statistic standard deviate = 2.0908, p-value = 0.03654
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.120686124     -0.014465249      0.004178391
chi.poly@data$chi.ols.res<-resid(chi.ols) #residuals ols


spplot(chi.poly,
       "chi.ols.res",
       at=seq(min(chi.poly@data$chi.ols.res,na.rm=TRUE),
              max(chi.poly@data$chi.ols.res,na.rm=TRUE),
              length=12),
       col.regions=rev(brewer.pal(11,"RdBu")))

Including Plots

You can also embed plots, for example: